#This script executes event-study regressions to estimate cohort-specific 
#associations between first-generation prenatal exposure to nonattainment 
#and different measures of TSP exposure for first- and second-generation individuals.

rm(list=ls())
gc()

library(reldist)
library(ineq)
library(stargazer)
library(data.table)
library(dplyr)
library(dtplyr)
library(ggplot2)
library(stringdist)
library(lfe)
library(haven)
library(readr)
library(broom)
library(tidylog)
library(rdrobust)
library(HonestDiD)
library(did)
library(here)
library(fixest)

# Set the working directory
setwd("/projects/opp_env/caa1970/")

# Source external R scripts for custom functions
source("Code/Child_Outcomes/child_outcomes_functions.R")
source("Code/rounding.r")

# Load the second-generation analysis dataset
load("Data/second_gen_college_analysis_data_jpemicro.rda")

regdata_rr <- regdata %>% filter(year%in%1969:1975, !is.na(all_fullterm), !is.na(parent_race), !is.na(PI_PC))

prenatalTSP_nocontrols <- feols(all_fullterm~ i(year, nonattainment_CAA70, ref=1971) | #Controls
                                  state_year+county_factor+year_factor+month+
                                  survey_year_by_year, cluster = "county_factor", data=regdata_rr )
betahat <- summary(prenatalTSP_nocontrols)$coefficients
sigma <- summary(prenatalTSP_nocontrols)$cov.scaled

delta_RM_results <- HonestDiD::createSensitivityResults_relativeMagnitudes(betahat = betahat,
                                                                           sigma = sigma,
                                                                           numPrePeriods = 2,
                                                                           numPostPeriods = 4,
                                                                           Mbarvec = seq(0.1,1,by=0.1),
                                                                           alpha = 0.05
)

write.csv(delta_RM_results, file="JPE_Micro/Output/Figure 1/delta_RM_results.csv", row.names=F)
